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Abstract 



The vast majority of recent advances in the field of numerical radiative transfer 
relies on approximate operator methods better known in astrophysics as Accel- 
, erated Lambda-Iteration (ALI). A superior class of iterative schemes, in term of 

^ ' rates of convergence, such as Gauss-Seidel and Successive Overrelaxation methods 

■ were therefore quite naturally introduced in the field of radiative transfer by Tru- 

i jillo Bueno and Fabiani Bendicho [1]; it was thoroughly described for the non-LTE 

] two-level atom case. We describe hereafter in details how such methods can be gen- 

eralized when dealing with non-LTE unpolarised radiation transfer with multilevel 
. atomic models, in monodimensional geometry. 

^_ 
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1 Introduction 



Our ability to deal with complex radiative transfer problems has considerably 
improved during the last ten years. It is indeed still worth an effort since im- 
portant issues of astrophysical interest may require to address a number of 
cases ranging from multi-dimensional problems in various geometries (Auer 
and Paletou [2]; Auer, Fabiani Bendicho and Trujillo Bueno [3] ; Fabiani Ben- 
dicho, Trujillo Bueno and Auer [4]; van Noort, Hubeny and Lanz [5]; Gout- 
tebroze [6]) to polarised radiation transfer involving complex atomic models 
(Trujillo Bueno and Manso Sainz [7]; Trujillo Bueno [8]), for instance. 

Most of the recent developments in the field of numerical radiation transfer are 
based on the Approximate (or Accelerated) Lambda-Iteration (ALI) scheme 
which has been recently reviewed in [9]. ALI methods are based on operator 
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splitting and, in a seminal article, Olson, Auer and Buchler [10] demonstrated 
the merits of adopting an approximate operator which is nothing more than 
the exact diagonal of the full operator A allowing the determination of the 
radiation field from a known source function (see Mihalas [11]; Rutten [12]). 
More generally speaking such an efficient iterative scheme is better known as 
the Jacobi's method in numerical analysis (Young [13]). 

However, even if superior - in term of higher rates of convergence - iterative 
schemes based on the Gauss-Seidel (GS) and the Successive Overrelaxation 
(SOR) methods have been introduced since in the field of numerical radiative 
transfer, they still deserved too little attention by the astrophysical commu- 
nity. GS/SOR methods have been described in details in a landmark article [1]; 
in particular, the practical description of how to implement GS/SOR itera- 
tions within a regular short characteristics formal solver is very meticulously 
written and therefore extremely useful. However, their description of the im- 
plementation of GS/SOR methods was restricted to the two-level atom case 
in monodimensional (ID) geometry. 

In a subsequent article [4], although mentions to GS/SOR iterative schemes 
generalized to the multilevel atom transfer problems can be found, their im- 
plementation is far from being made explicit and it is only described in very 
general terms; it is furthermore embedded in another (very efficient though) 
numerical strategy based on multi-grid methods [14,15]. Unfortunately, it ap- 
pears that the basic features of a GS/SOR-based formal solver in the frame of 
the multilevel atom radiative transfer problem remain, so far, to be explicited. 

The present article aims indeed at filling the gap by providing all the ele- 
ments required for a successful implementation of multilevel GS/SOR iter- 
ative schemes in ID geometry. We shall recall in §2 the basic principles of 
GS/SOR iterative schemes in the case of a two-level atom model. Then, in 
section §3 we shall describe step by step how GS/SOR can be implemented 
for the case of multi-level atom models, therefore extending the well-known 
multilevel- ALI method of Rybicki and Hummer [16]. Finally, we shall present 
in §4 some illustrative examples clearly demonstrating the performances of 
multilevel GS/SOR iterative schemes. 



2 Gauss-Seidel and SOR iterative schemes basics 

Gauss-Seidel and SOR iterative scheme are particularly well adapted to the 
short characteristics method (SC) for the formal solution of the radiative 
transfer equation [2,17,18]. 

Using SC in ID geometry indeed, the formal solution is obtained by sweeping 
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the grid say, first in direction —fl [fj, < 0) i.e. from the surface down to 
the bottom of the atmosphere, and then in the opposite, upward direction 
+fl (/i > 0) starting from the bottom of the atmosphere up to its surface 
though. The specific intensity Im is advanced step by step during each pass, 
partially integrated over angles and frequencies during the downward pass 
while, during the second (upward) pass, the mean intensity J can be fully 
computed, completing therefore the formal solution at each depth 



Jk = MSk] . (1) 



Except at the boundaries where the illumination conditions are known a priori, 
along each direction, the specific intensity at the inner grid points is advanced 
depth after depth, and computed along the short characteristics (allowing us 
to suppress any index related to the angle variation of the specific intensity 
hereafter) according to 



lo = 4e-^-" + ^^Su + ^oSo + "^dSd (2) 



where the first part of the right-hand side of this expression corresponds to the 
part transmitted from the "upwind" grid point u down to the current point o 
(see Fig. 1), and the three last term result from the analytic integration of 



T = y S{T)e-^dr (3) 



along the short characteristics going from u to o. Indeed, assuming that S is 
quadratic in the optical depth r, it is easy to show, using a simple Lagrange 
polynomial interpolation on the basis of three consecutive grid points u, o and 
d, that 



T^^^S^ + ^H^So + ^dSd, (4) 



3 



where the ^'s coefficients are defined^ as 



^ u 



W2 + WiATd 

At„(At„ + Ara) 



Wi{ATu - ATd) - W2 



d 



and where 



Wq 



W2 - WiATu 

Au{Atu + ATd) 



1 - e-'^^" 



Wi = Wo- Ar„e~^^" • (6) 



W2 



= 2wi - Ar^e"^^" 



In the two-level atom case, the non-LTE fine source function S is usually 
expressed as 

S{T) = {l-e)J{T)+eB{T), (7) 



where r is the optical depth, e is the collisional destruction probability mea- 
suring departures from LTE, B the Planck function and J is the usual mean 
intensity 

dQ 



where we omitted the optical depth dependence for simplicity. 

The Jacobi type iterative scheme introduced in numerical radiative transfer 
by [10] allows for the iterative determination of the source function on the 
basis of increments such that 



l-(l-e)A 



kk 



^ A different expression was used in [2], Eqs. (8) to (10); note also that there is a 
sign error in their coefficients di and d2 ■ But more important is to adopt the current 
expansion in terms of ^''s for the efficient implementation of GS/SOR. 
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Fig. 1. Three successive grid points along direction +Q (upward) are used, assuming 
parabolic interpolation of the source function, in order to compute the specific 
intensity at current position o along the short characteristics starting from u. 

at each depth Tk, where the superscript (old) denotes quantities already known 
from the previous iterative stage and where A^.^ is a scalar equal to the diag- 
onal element of the full operator A at such a depth in the atmosphere. 

Within a Gauss-Seidel iterative scheme the source function increments now 
turn themselves into 

/I \ rfold and new) , 7-, r((old) 

A niGS) _ (l-£)4 '+eBk-Sl 
1 - (1 - e)Akk 



where the somewhat enigmatic quantity means that at the spa- 

tial point k the mean intensity has to be calculated via a formal solution of 

the transfer equation using the "new" source function values Sj^"^^^ already 
obtained at points j = 1, 2, . . . , {k — 1), and the "old" source function values 
Sf'^'^ at points j = k, {k + 1), . . . 

Achieving this was explained in every detail in [1] and we shall not repeat 
here the algorithmics of the process. It is however important to repeat that 
such a scheme is better implemented in the frame of the short characteristics 
approach since it allows to avoid any explicit matrix inversion. 

Finally, going beyond Gauss-Seidel, SOR iterations can be implemented where 

A5f = cA^^^) , (11) 
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uj being, in the case of overrelaxation, a parameter such that 1 < cu < 2, with 
an optimal value ~ 1.5 for the two-level atom case [1]. 

However, in the non-LTE multilevel atom case, the radiation transfer problem 
is somewhat more complicated since it is now required to solve self- consistently 
for the coupled set of A^trans radiative transfer equations, for all radiatively 
allowed transitions, together with a set of A^icvcis equations of the statistical 
equilibrium (ESE) at each optical depth giving the populations numbers for 
all atomic levels. 



3 The multilevel atom case 

Hereafter, we shall describe in details how the GS/SOR numerical method 
can be implemented for the case of multilevel atom models in ID geometry. 
For simplicity wc shall consider the case of non-overlapping lines with no 
background continuum. 

3.1 The MALI method: an overview 

ALI methods have been first generalized to the multilevel atom case by Ry- 
bicky and Hummer [16,19]. Their MALI method consists in the preconditioning 
of the equations of statistical equilibrium (ESE) using an approximate oper- 
ator A*j and a J^f mean intensity defined hereafter, and computed, for all 
allowed transitions u ^ I, from quantitites known from the previous iterative 
step. 

Generally speaking, the line source function for a radiative transition between 
the two bound atomic levels u (i.e., the upper energy level) and / (respectively, 
the lower one) can be written in terms of the density of populations nu,i, the 
Einstein [20] coefficients for spontaneous emission A^i, absorption and 
stimulated emission Bui, and the line absorption (f)^, and emission if^i, profiles 
such as 



Hereafter we shall however restrict our study to complete redistribution in 
frequency for which absorption and emission profiles are identical i.e. 0^ = ijji, 
and therefore, the line source function remains independent of the frequency. 
It is indeed well-known that such a assumption remains valid for most of the 
spectral lines; however, partial redistribution in frequencies effects should be 




(12) 
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considered, at the cost of more computational work, for a proper description 
of resonance line formation in a diluted medium (e.g., [22] and references 
therein). 

Furthermore, leaving out bound-free transitions, at each depth in the atmo- 
sphere, the ESE are usually expressed as a set of Nieyeis elementary equations 
such that 



j<i 

~y^. ['T'j^ji ~ {''^iBij — njBji)Jij] 

j>i 

+J2i^iCij-njCji) = 0, (13) 
j 

where the Cij are collisional excitation/dccxcitation rates. Now, introducing 
the approximate operator into the ESE via the ALI approximation 

I^a ^ Kn[S] + (Kn - A:n)[^^°"^] , (14) 



where S^°^^^ represents the source function known from previous iteration and, 
defining 

(15) 



where J^^ depends only on known quantities, we can establish ^ after [16] the 
following set of preconditioned equations 



j<i 
j>i 

+J2i^iCij-n,C,i) = 0, (16) 
j 

where 

A-= J^J<P.Kndi^. (17) 

^ Inserting Eq. (12) into Eq. (14) and forming the net radiative rate {rijBji — 
niBij)Jij in terms of A*- and J?f. 
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Heinzel [21] and Paletou [22] showed how to take care self-consistently of 
the ionization equihbrium within the MALI approach by adding a Newton- 
Raphson scheme to it, for the determination of the electron density from a set 
of non-hnear ESE taking into account all kind of bound-free transitions; the 
same can be done within the new GS/SOR multilevel numerical schemes. 



3.2 Expliciting GS/SOR with multilevel atoms 

Assume that one has already swept the grid once, say from the illumination-free 

surface of the atmosphere at A; = 0, down to the bottom boundary at A; = ND 
along direction —Q. By analogy with the GS/SOR numerical strategy for the 
two-level atom case, in the multilevel atom case we are now going to update 
all population numbers at successive depths k — ND, 1 while sweeping back 
the grid along the opposite, upward direction 

Now the population update will be made depth after depth, from the bottom, 
up to the surface of the atmosphere by inverting the MALI preconditioned set 
of ESE given in Eqs. (16) before passing to the next depth point. It is a quite 
straightforward task at the lower boundary surface since the incident radiation 
field is known a priori from the (given) external conditions of illumination. 

The situation is however a bit more tricky at the inner grid points. Once the 
populations at depth {k + 1) have been updated, we shall advance along direc- 
tion +Q to the next grid point at depth k. But having changed {nj'^'^^}(^k+i) to 
{n^°'^'^^}(fc+i) means that the local absorption coefficients X(fc+i) and the source 
functions for all allowed transitions have to be changed accordingly, as well as 
the upwind optical depth Ar^l]^^^ along the path from depth (k + 1) to depth 

k. As a consequence, the three coefficients [d = {k — 1), k, {k + 1)] used 
for the evaluation of the specific intensity along the short-characteristics are 
also affected by the local population change and thus need to be updated. But 
since At^I^^-^-^ = Ar^^^^j^^ evaluations of the specific intensities made during the 
first downward pass must also be corrected accordingly, for consistency. 

So the major point in such implementation is to carefully "propagate" the 
effects of the local population update. Let us, for instance, consider the specific 
intensity evaluation at any inner grid point k in, respectively, direction —fl 
for which we shall use (j) superscripts, and +^ (resp. (t) superscripts); we 
can therefore write, following Eq. (2) 
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1 e-M^ — • ^ ' ^ ' ^ ' ' ' ^ ' ' ' ^ 

100 200 300 400 500 600 700 

iteration number 



Fig. 2. Rates of convergence for the MALI, Gauss-Seidel and SOR multilevel iter- 
ative processes, respectively. A spatial grid of 20 depth points per decade together 
with a 8 Gauss-Legendre angular quadrature and constant Doppler profiles were 
used. The atomic model is a 3-level H l model atom taken from [23] (see also Table 
2). 

while, at the same depth but in the opposite direction, we have 

4" = + *m,v, + + *<",S._, . (19) 



Since the update of the population numbers at depth (A; + 1) have just been 
done by inverting the system of Eqs. (16), this generates for each allowed 
transition changes in the absorption coefficients at line center, now becoming 



(new) 
^(fe+1) 



n 



(new) 
7 



lu 



ul 



(20) 



and in the line source functions, assuming complete redistribution in frequency, 
turning into 

Ci(new) ' ''-u ^ul /i^^ \ 



Therefore, one has first to correct Ij^ ' for consistency since previous changes 
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Fig. 3. Source functions vs. optical depth for the three transitions of our H-model. 
Sohd Hnes are our GS/SOR multilevel results and dashed lines correspond to 
low-resolution results previously published in [23] and [26]. 

lead, according to Eq. (18), to 



.a),(new) _ .a), (old) -Ar^i^' ^(i), (new) ^(old 



Did) 



_|_ ^(i). (new) ^(old) ^(i), (new) ^(new) 



(22) 



which requires that ij^^l'i"^^^ have been saved in memory during the first down- 
ward pass. This step is equivalent to the computation of the AJj^ correction 
pointed out in [1], Eq. (39). 

Then, ij^^^ can be computed in a straightforward manner via Eq. (19), which 
makes it possible to compute A*j and J^? for all allowed transitions, then to in- 
ject these quantities in the preconditioned ESE and finally to compute/update 
locally {nj'^^^^lfe while in the middle of the upward pass. 

However, before advancing to next depth point (k—l), it is important to con- 
sider first the various changes induced by the level population changes respec- 
tively on the source functions Sk, on the optical depth increments Ar^^^^^p^^ 

and AtI^I^_^[^^'"^ as well as, finally, on the three short characteristics coefficients 

^J) [d= {k-1), k, {k + 1)]. Because of that, we end up all updates at depth 
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k during the upward pass by the computation of 

(23) 

_|_ ^(T), (new) ^(new) _|_ ^(T), (new) ^(old) 
k k k — 1 k — 1 

This last stage is analogous to the correction described by Eq. (40) in [1] . 

Finally, a multilevel SOR iterative scheme is built when, at each depth k, all 
the populations of the excited levels are updated according to 

^(new)^^(old)^^^^(GS)^ (24) 

where a; is a parameter of the order of 1.5 as discussed below. 

3.3 Numerical recipes 

Several "recipes" should be followed when implementing a GS/SOR solver for 
multilevel atoms. The first one is to order properly the various loops. From 
outer to inner loops one may find indeed: (1) the directions (±^^) along which 
the slab will be swept, (2) the number of allowed radiative transitions, (3) the 
direction cosines (i.e., the usual //'s) and, finally (4) the frequencies. 

Upwind and downwind corrections, as described in Eqs. (22) and (23), require 
some bookkeeping of variables such as all the ij^^^' after the downwind pass 
— Q for the further computation of the mean intensity entering the precondi- 
tioned ESE i.e., Eqs. (16). 

Table 1 

Computation time for the 3-level H model of [23] obtained with an Intel Pentium-4 
clocking @ 3 GHz as a function of the number of depth points per decade Nr', the 
stopping criterion was Rc = 10^^° (we also indicate the # of iterations required 
respectively) . 



Nr 


MALI 


GS 


SOR 


5 


7.644s (160) 


4.595s (78) 


2.325s (39) 


10 


27.852s (326) 


17.148s (161) 


6.096s (55) 


15 


59.768s (483) 


38.997s (239) 


10.203s (63) 


20 


lm43.373s (634) 


lm05.879s (315) 


17.213s (81) 


25 


2m37.637s (780) 


lm40.630s (387) 


26.040s (99) 
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400 



200 

iteration number 

Fig. 4. History of the convergence error (dashed hnes) and of the relative change 
Rc (fuh hnes) for, respectively, the GS and the SOR multilevel iterative schemes. 
Although a small value of Rc does not necessarily imply a small enough value of Cg as 
demonstrated by the GS curves, a SOR iterative process including a self-consistent 
evaluation of uj leads to a scheme such that Rc — Ce- 

Finally, it is important to realize that the additional computing time necessary 
for all the extra-computations required by one GS/SOR iterative step is very 
small as compared to the huge saving on the whole convergence process with 
respect to MALI (see Table 1). 



4 Illustrative examples and discussion 



We have adopted the standard benchmark models for multilevel atom prob- 
lems proposed by Avrett [23] and Avrett and Loeser [24] whose main param- 
eters have also been summarized in Table 2. 

The respective rates of convergence for the MALI, Gauss-Seidel (GS) and SOR 
multilevel iterative processes are displayed in Fig. 2 where we plotted the max- 
imum relative change on the level populations from an iteration to another, 
Rci for the schematic 3-level Hi model. A spatial grid of 20 depth points per 
decade together with a 8 Gauss-Legendre angular quadrature and constant 
Doppler profiles were used. It clearly demonstrates how superior to MALI the 
multilevel SOR iterative scheme can be for such refined grids; and it is worth 
pointing out here again the importance of grid refinement on the accuracy of 
ALI-class methods, as demonstrated by Chevallier et al. [25]. Indeed, in such 
factor of ~ 6 in computing time can be saved, when one iterates 
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Table 2 

Input parameters for the Hi and Can multilevel benchmark models of Avrett [23]. 
Statistical weights are gi = 2, 52 = 8, and 33 = 18 for Hi, gi = 2, g2 = 4, 53 = 6, 
g4 = 2, g5 = 4 for Call; the temperature of the atmosphere is 5000 K. 



element 


u 


/ 


Aui 






H 


2 


1 


4.68 X 10^ 


10^ 


2.47 X 10^^ 


H 


3 


1 


5.54 X 10'^ 


10^ 


2.93 X 10^5 


H 


3 


2 


4.39 X lO'^ 


10-'^ 




Ca 


2 


1 


forbidden 


8.2 X 10^ 


4.10 X 10^"^ 


Ca 


3 


1 


forbidden 


8.2 X 10^ 


4.12 X lO^'' 


Ca 


4 


1 


1.4 X 10^ 


5.1 X 10^ 


7.56 X 10^^ 


Ca 


5 


1 


1.4 X 10 


0.1 X lU 


7.00 X lU 


Ca 


3 


2 


forbidden 


10^ 




Ca 


4 


2 


7.8 X 10^ 


1.6 X 10^ 




Ca 


5 


2 


8.1 X 10^ 


1.6 X 10^ 




Ca 


4 


3 


forbidden 


10'^ 




Ca 


5 


3 


7.2 X 10^ 


1.4 X 10^ 




Ca 


5 


4 


forbidden 


4.8 X 10^ 





the population numbers down to numerical noise, for instance. We report a 
similar behaviour of the iterative processes when dealing with Avrett's Call 
ion model. 

In Fig. 3, we display (solid hnes) the source function vs. line optical depth for 
the 3 transitions allowed by our schematic H-model. Dashed lines correspond 
to the old benchmark results of [23]; differences are coming not only from 
very different numerical schemes but also from large differences in the angular 
and frequency quadratures. Our results are also in good agreement with those 
of [26] and may serve, given the high-level of grid refinement we adopted, as 
new benchmark results for multilevel atom cases. 

The main potential drawback of SOR methods is that it relies on the choice of 
a relaxation parameter u whose optimal value is a priori unknown. However, it 
was proposed [1] a quite robust numerical procedure in order to estimate self- 
consistently a close-to-optimal cu after having run a few "pure GS" iterations. 
We followed their recommendations and, indeed, found a posteriori values 
close to "optimal" ones deduced from experimental runs using a prescribed u>. 

Finally, in Fig. 4, we plotted the history of the convergence error Cg defined 
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as 



Ce = max 



( 



n{itr) — n(oo) 



n{oo) 



) 



(25) 



following [3], where itr is the iteration number and n{oo) is the fully converged 
solution, and of the relative change Rc for, respectively, the GS and the SOR 
multilevel iterative schemes. It is important to note again that reaching a small 
value of Rc, which is indeed the most direct control parameter of the iterative 
process, does not necessarily imply a small enough value of Ce to guarantee 
convergence; and this is shown by the GS curves. However, an optimal SOR 
iterative process including the self-consistent evaluation of u as proposed in [1], 
leads to a better-controlled process since Rc is just slightly lower than Ce, as 
shown in Fig. 4. 

The issue of having a reliable stopping criterion is of course critical to any 
iterative method and, for the specific case of numerical radiative transfer, 
Auer et al. [3] and Fabiani Bendicho et al. [4] addressed it successfully by 
adopting multi-grid methods. We finally refer the reader to this later work 
where it is demonstrated how the combination of GS/SOR iterative schemes 
together with multi-grid techniques lead to extremely powerful techniques for 
the solution of complex radiative transfer problems. 
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